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Abstract 

We describe a new microscopic approach for analyzing interacting electron 
systems with local moments or, in principle, any local order parameter. We 
specialize attention to the doped Mott insulator phase of the Hubbard model, 
where standard weak-coupling perturbation methods fail. A rotationally in- 
variant Stratonovich-Hubbard field is introduced to decouple the static spin 
components of the interaction at each site. Charge degrees of freedom are then 
treated by "slave" Hartree-Fock in the presence of a spatially varying random 
spin field. The effective action reduces to a classical Heisenberg model at half 
filling, insuring that the system has (i) finite-range order at T > with an expo- 
nentially diverging correlation length and (ii) a one-electron Mott-Hubbard gap 
in the presence of disorder. Away from half-filling properties are determined 
by strongly non-Gaussian fluctuations in the amplitude and orientation of the 
local spin fields. The saddle point equations of the theory at zero temperature 
reduce to inhomogeneous Hartree-Fock, so that disordering of domain walls 
at finite temperature is in principle included. We present preliminary small 
system results for the intermediate and large interaction regimes obtained by 
Monte Carlo simulation of random spin field configurations. 
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1 Introduction 

In this paper we describe a new approach to perturbation theory for systems with 
local moments or, in principle, any local order parameter. We begin by summarizing 
the present status of weak-coupling perturbation theories for such systems, with the 
Hubbard model as a primary example. 

Over the course of the last decade weak-coupling perturbation theories have been 
used extensively to investigate the physics of the Hubbard model and its extensions. 
The simplest perturbative approach to the Hubbard model is just Hartree-Fock theory 
for one-body correlation functions. The conserving approximation implied for two- 
body correlations is then the random phase approximation, or RPA. In the large-?/ 
limit at half-filling the static spin response is unstable within Hartree-Fock, and the 
system exhibits a broken-symmetry ground state with infinite-range antiferromag- 
netic order [0. Furthermore a gap of order U appears in the one-body density of 
states for the ordered phase. If the Hartree-Fock equations are extended to finite 
temperature, long-range antiferromagnetic order onsets at a temperature of order U 
in the large-?/ limit. The zero-temperature behavior of the Hartree-Fock approxima- 
tion away from half-filling has also been investigated, and incommensurate states, as 
well as vertical and diagonal domain walls, have been studied . However, the appro- 
priate treatment for finite temperatures has remained unclear, because the energetic 
competition of the various states is delicate. 

In fact the Hartree-Fock approximation at finite temperatures provides only a 
crude description of the physics in a two- (or one-) dimensional system. In accordance 
with the Mermin- Wagner Theorem [[| , infinite-range magnetic order is forbidden for 
T > in two dimensions. The magnetic correlation length £(T) instead diverges 
exponentially in 1/T for T — > at half-filling, and the quantitative behavior of 
£, as well as the zero-temperature ordered moment, is renormalized by quantum 
fluctuations. Nevertheless the half-filled model exhibits a gap (or pseudogap, since 
we do not mean to imply a rigorously zero density of states) at finite temperature in 
the spin-disordered phase. 

With the intention of modeling strong spin fluctuations in a disordered system, the 
present authors proposed the fluctuation exchange, or FLEX, approximation j|, |, || 
as a plausible starting point. The self-consistent treatment of spin fluctuations in 
the one-body self-energy suppresses the magnetic ordering transition to zero tem- 
perature. Nevertheless the new approximation is itself flawed: It does not produce a 
sharp gap at half-filling and instead smears out the density of states, leaving a finite 
weight at zero energy for all temperatures |5|, ^|, 0. Furthermore the spin response 
function diverges with decreasing temperature in a way which only roughly mirrors 
the exact behavior of the model. The failure at half-filling calls into question results 
obtained for large U at small doping, and in the end the FLEX approximation fails 
to describe the pseudogap regime. 

The inadequacy of the FLEX approximation has motivated various attempts to 
improve the theory. Kampf and Schrieffer have calculated the electron spectral density 
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near an antiferromagnetic transition to lowest order in an effective spin-fluctuation 
interaction, and their results exhibit a pseudogap when the spin correlation length is 
large M. Vilk and Tremblay || have introduced a factorization approximation for 
the electron four-point function which is constrained to be exact when all space-time 
coordinates coincide. Others |TD| have argued that perturbative calculations based 



on an unrenormalized electron propagator are justified by a cancellation between 
vertex corrections and self-energies, leading to the appearance of a pseudogap in 
some parameter regimes. 

In addition efforts have been made to extend the FLEX approximation within 
a scheme based on infinite-order (or "parquet") vertex renormalization O, |T^|, |T3|| . 
Quantitative studies of parquet renormalization have had the goal of establishing the 
validity or inadequacy of perturbation theory resummations for systems with a gap, 
but without infinite-range order. Why would one assume that such resummations 
have any chance of succeeding? For at least one nontrivial model, the Anderson impu- 
rity in a metallic host, it has been known since the mid-1980's that zero-temperature 
perturbation theory for thermodynamic properties has an infinite radius of conver- 
gence in U |L4|. This fact, as well as some suggestive results in the low-density limit, 
initially encouraged hope that an infinite-order parquet summation for the doped 
Mott insulator might succeed if the technical difficulties impeding its implementation 
could be overcome. 

In fact this hope is not justified. Full numerical solutions of the minimal parquet 
equations for the two-dimensional Hubbard model in the large— U limit [TjJ yield 
the following results: (i) No gap appears in the converged solution at half filling, 
(ii) The spectral density at low energies is actually larger than that in a FLEX 
calculation. We interpret these results to mean that conventional perturbation theory 



bj] fails to describe the physics of a Mott-Hubbard gap or pseudogap in a spin- 



disordered system. Furthermore the parquet results, which consistently incorporate 
vertex renormalization to infinite order, bring into question whether the basic problem 
can be circumvented by the constrained factorization approximation or a cancellation 
between parts of perturbative diagrams. 

Rather we believe the failure of perturbation theory around a uniform (i.e., Fermi 
liquid) starting point is due to a qualitative change in the saddle-point structure 
of the action for the Hubbard model at temperatures of order U in the \arge-U 
limit. This change is induced by the appearance of local moments (or, more generally 
speaking, a local order parameter field) in the theory. Conventional perturbation 
theory becomes unreliable below a temperature scale which depends on both the 
interaction U and the chemical potential //. Below this scale one needs a revised 
perturbation theory which incorporates "anomalous" Feynman graphs, but does so 
without inducing infinite-range order. These seemingly contradictory objectives may 
be accomplished by introducing a spatially varying static field with the full symmetry 
of the associated classical order parameter. One way to introduce such a field is via a 
Stratonovich-Hubbard (S-H) transformation [17 , 18 of the interaction. In particular 
we employ below a spin-rotation-invariant S-H field to decouple the static (m = 0) 
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component of the Hubbard interaction. About thirty years ago Wang et al. [19 
showed for the magnetic alloy problem that the effective static free energy induced 
by a time-varying S-H field provides a physical picture of the crossover from band- 



like to local-moment behavior. Evenson et al. |2(| analyzed the coupling between 



isolated moments which arises within this approach and noted that the simplest 
scheme generates an Ising, rather than Heisenberg, interaction. Below we show that 
when a vector S-H field is introduced, the effective action reduces to a classical 
Heisenberg form at half-filling. This means no infinite-range order at T > 0, an 
exponentially diverging correlation length £(T), and a Mott-Hubbard gap in the spin- 
disordered phase. Furthermore the approximation scheme should provide a systematic 
way to study doping of the Mott insulator. In the doped case both orientation and 
amplitude fluctuations of the field become important. In particular the effective action 



for the static field is strongly non-Gaussian | 21[| . Since the action arises directly from 



integrating out the fermion degrees of freedom, its form reflects both the underlying 
spin and charge structure of the system. In this way domain wall stripe fluctuations 
can occur naturally if they are present in the model. 

We describe this new approximation scheme for the Hubbard model in Section 
2, then derive expressions for specific electron correlation functions in Section 3. In 
Section 4 we provide a brief description of our Monte Carlo sampling scheme for 
integration over the fluctuating local fields. In Section 5 we present some illustrative 
results obtained on a work station for a 4 x 4 system. We outline in Section 6 the 
form of generalized approximations which extend the present approach in the same 
sense that FLEX extends Hartree-Fock. Finally we close in Section 7 with some 
speculations on possible implications for high-temperature superconductivity and 
with an outline of the next steps in this work. 



2 Approximation scheme 

In this section we show how to decouple the the static spin interaction in the 
repulsive Hubbard model by introducing a spin-rotation-invariant S-H field. This 
field may be interpreted as a classical Heisenberg order parameter which generalizes 
the saddle-point in conventional temperature-dependent Hartree-Fock calculations. 

It is conventional to decouple the interaction term by introducing independent 
fields at each point in imaginary time. When this is done, some freedom remains in 
choosing the form of the interaction, due to the equal-time fermion anticommutation 
relations. One may write equivalently 

Ufkfa = \Uinl-ml) (1) 
= \Uirii-ml) (2) 
= \Un i -\U{ml + rn%) , (3) 

where 

% = n iT + n a 
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There is in fact considerably more freedom allowed, since the interaction is rotation- 
ally invariant: The z-component of the moment operator in Eqns. [I] and || may be 
replaced by the x- or y-component, and likewise different pairwise combinations of 
moment operators may be substituted in Eqn. |3|. One is free to consider in addition 
arbitrary linear combinations of these expressions, weighted properly to reproduce 
the original interaction. 

This was realized by several authors around 1970 during early studies of the 
single-impurity Anderson model and the Hubbard model. In particular Wang et 
al. employed Eqn. |^ |TP[ and Hamann employed Eqn. [I] [25| as the starting points for 



equal-time Stratonovich-Hubbard decoupling schemes. Note that both these expres- 
sions imply an Ising-like S-H field governing the spin degrees of freedom, while the 
third expression implies an xy-like field. All three expressions result in actions which 
break the spin-rotational symmetry of the theory, and the symmetry is only restored 
by the integration over field configurations. As long as all configurations are included 
there is no difficulty, but problems arise if approximations are introduced for any 
system beyond a single impurity. This was realized by Evenson et al. [BO] in large-t/ 
studies of the two-site Hubbard model within a static approximation: The effective 
spin-spin coupling assumes an Ising, rather than Heisenberg, form when Eqn. BJ is 
used as a starting point. It would clearly be desirable to introduce a symmetrical 
Heisenberg decoupling, but this is difficult for the following reason: No matter how 
the spin quantization axis is chosen, the three equivalent components of the spin- 
spin interaction are divided up between the two particle-hole channels (longitudinal 
and transverse) of the electron vertex, which "overlap" completely at equal time. A 
Heisenberg decoupling can be introduced artificially if, for example, the interaction 
is written as an equal-weight linear combination of Eqn. |I] for the three choices of 
moment operator, i.e., 

Un^riii = ~Uf% — ^Umi-mi . (5) 

This procedure introduces spurious factors of | in the saddle-point equations when 
a direction is chosen for the broken-symmetry spin axis, however, and so offers an 
unsatisfactory solution. 

We may restate the problem as follows: A conventional equal-time S-H transfor- 
mation, followed by a static approximation within the resulting theory, violates the 
crossing symmetry of the Hubbard interaction vertex and prevents a spin-rotation- 
invariant treatment of fluctuations. To insure an equivalent static treatment of all 
three spin components, we propose the introduction of a zero-frequency, rather than 
equal-time, S-H field. Such a decoupling scheme is possible because the static compo- 
nents of the direct and exchange channels in the Hubbard model have essentially zero 
overlap at low temperatures. This allows the simultaneous treatment of interactions 
in both crossed channels and leads to a unique prescription for the decoupling scheme 
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independent of the expression employed for the interaction in Eqns. |T}|^. The same 
procedure can be applied in more general situations when it is necessary to decouple 
multiple interactions in crossed channels. 

It is possible to carry out this transformation directly within an operator-based ap- 
proach, but the manipulations are more transparent with anticommuting c-numbers, 
and we follow the latter approach. We consider a discretized representation of the 
partition function 

Z = Tie-?" , (6) 

where 

H = ~tJ2 c]+6,* C ia + U J2 C l C iAl C H ~ ^H C t C ia- (7) 

iScr i ia 

Let 

r e = £At 

At = P/L, (8) 

with L the number of time slices. We write the anticommuting c-numbers for time 
T£, site i, and spin a as c l f, cf . It is convenient to use the frequency-transformed 
variables c l °, c™, where 

L/2-1 

-ia _ e i^nT£ -^ia 

n=-L/2 

cf = ]T e-**"<% (9) 

n=-L/2 

with 

u n = (2n + 1)ttT (10) 
for integer n. In order to keep track of discretization effects, we define 



which approaches iu n for Ar — > 0. 
The discretized fermion action is 



ina i8na 

+ flU ^ c n+m c -n-l c -n'-l C n'+m ' (^) 
imnn' 

The interaction term is represented in Figure 1. (Note that the integer label — n — 1 
corresponds to frequency ui- n -i = —w n .) Equivalently one may write 

S = Sq + Si nt , (13) 
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with 



0£# 



-iu£l + (H - pi) 



c JrT 



(14) 



i] no- 



where the matrix in brackets has spatial indices, and 



"-i.i 



0)ij 



—t, i and j near-neighbors 
0, else. 



(15) 



Now note that the interaction can be rewritten in two equivalent ways to emphasize 
the particle-hole channels (see Figures 2a-b): 



Si 



hit 



fiU £ C n+m C n C n' C n' 



pu £ 



(16) 



At this stage we restrict attention to terms with m = 0, i.e., the static component 
of the particle-hole interactions. For m = the two sums in Eqn. |l^ are essentially 
independent: The only terms which appear in both sums are restricted to n' = n. 
Therefore the static particle-hole interaction may be written 

= -#/£[£##] [£#£,] + ^£[£4 T 4 T ][£^4] (it) 



where we have dropped from the double sums an n! = n correction term, which 



becomes negligible at low temperatures [23 



Now all terms in Eqn. [17] may be reduced to sums of perfect squares in preparation 
for the introduction of S-H fields. First note that 



\ £ < a K 4 irr 



, , c ia 



\ £ C ipx - w 



, , c ia 

y ) aa 



Using these identifications one has 



[ E C n C n ] [£ C n' C n> 



(<o) 2 + Ko) 2 



(19) 



where 



in 



.r0 



£ "" ! 

naa 



yO 



naa' 



^n \yx jaa' ^r, 

-Jta t \ it 
c n \yy)am' c n 



(20) 



s 



Likewise, 



n 



where 



??. n 



m 



zf) 



IK-<o) 



( c n c n "I" C n C n ) 
X # 



(21) 



z J act' ^ n 



In this case one identifies 

[ X c n c n ] [ X/ S n' c n' 

Thus 



nl)' 2 - K ) 2 



m=0 
int 



m ■ m 



i \2 



(22) 



(23) 



(24) 



The repulsive charge interaction is retained for treatment by perturbation theory 
(i.e., saddle-point plus fluctuations). The attractive spin interaction in Eqn. |24] may 
be decoupled using an independent static S-H field at each site i: 



cxp 



-s. 



m=0 
spin 



cxp 



m ■ m 



J n(^) 3/2 c^;c ex P [-/3£f-f + pvuy:c 



m. 



(25) 



The S-H fields £ * may be interpreted as classical magnetic fields which couple to 
the local electron spin polarization. The Hubbard action may be rewritten as 



s = /3£ lei 2 + So 

i 

where 

So = P X Z% n { -^n^aa'Sij + (Ho)ij5 aa > 



nm=0 _|_ qm^O 
•^charge "l - '-'int ' 



ijnaa' 



and 



qm=0 
charge 



n, 



j \2 



(26) 

(27) 
(28) 



The interaction contributions to the action for m/fl take precisely the form they 
would in a conventional perturbative treatment, except that all vertices related to 
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the zero-frequency spin vertex by crossing are absent. Note that because we have 
decoupled only the m = spin components, we have left open the possibility of higher 
order perturbative expansions about a static approximation. 

We consider now the static charge interaction S^^ ge in Eqn. This term could 
be decoupled in parallel with S"™^ by introducing a real-valued S-H field p l . Since 
the charge interaction is repulsive, however, the fermion bilinear coupling to p % would 
be pure imaginary and unsuitable for Monte Carlo simulation. It is straightforward to 
show that a saddle-point approximation for the p % fields in the presence of a fixed £ 1 
configuration yields the same result as a spatially non-uniform Hartree-Fock. For this 
reason we do not introduce p l fields, but reserve S^rge f° r treatment by perturbation 
theory (i.e., Hartree-Fock at lowest order). 

Within such a Hartree-Fock approximation for fixed £\ 

ssSb - I^EWtW - i/^ £<**>?. (29) 

ina i 

where (ri 1 )^ is the mean occupancy at site i. It is essential to retain the "constant" 
term, since it depends on the local field configuration {£*}. From this point we use 
the notation (A)^ to denote the expectation value of the fermion observable A in the 
presence of a fixed field configuration. 

The static ^-dependent action assumes the form 

S(0 = S (0 + Sf(e, c, c) , (30) 

where 

MO = /?E lef - J/h/EM? 

i i 

S{(£,c,c) = j3 E c^^-iuj^dijdaaf + [(Ho);j(W - VUC ■ a aa >5ij 

+ (i[/(n^-//)(%(W]}c' ■ ( 31 ) 



ijnaa' 



The fermion contribution to the action may now be integrated out exactly, since it 
takes a one-body form (albeit a form non-diagonal in spin): 

f DcDcexp[-S { } = J] [l + , (32) 

J r 

where the e r are eigenvalues of the one-body Hamiltonian 



H(0 = -*E C ts,aC ta ~ V^E C-ctataa'C^ + E [l^Vf 



Thus the final expression for the effective action after the fermion integration is 

S eS (0 = 3„(0 - E log [1 + e-*-®] • (34) 
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Note that for fixed \x and a Hartree-Fock iteration is required to adjust the 
non-uniform charge field (n l ) to satisfy the condition 

<**> = E/WEiWP 

r cr 

/(^) = 1^7- ( 35 ) 

In this equation |r) is the one-body eigenvector corresponding to eigenvalue e r (£); 
and the \ia) constitute a complete set of basis vectors corresponding to spin S z = a 
on site i. Note that for an A^es x A" S ites real-space lattice, the size of the eigenvector 
space is 2A sites . 

The form of the charge-interaction Hartree-Fock employed here provides a tem- 
plate for the family of Baym-Kadanoff (or FLEX Q) generalizations of the 
present approach discussed in Section 6. The charge field is here viewed as a nonsin- 
gular "slave" to the spin field: It is fixed to its saddle-point, which changes when the 
£ configuration changes. 

No further approximation to the action S e s(0 is now possible without explicitly 
breaking spin rotational and translational invariance. In conventional Hartree-Fock 
treatments it is assumed that the static field assumes a rigid configuration. At 
half-filling the mean-field ansatz is just 

f = ge^* Q = (vr, tt) 
(n*) = 1 . (36) 

A saddle-point evaluation of S e s within this restricted field space reproduces the con- 
ventional Hartree-Fock equations for a commensurate antiferromagnet. Furthermore 
inhomogeneous Hartree-Fock solutions for (n) 7^ 1 follow by assuming a spin- 
space form (collinear, spiral, etc.) for the order parameter and solving the resulting 
saddle-point equations. 

We propose instead to evaluate electron correlation functions using S e s as the 
action in a classical field Monte Carlo. The Monte Carlo is classical since all finite- 
frequency contributions to the action have been set aside for treatment by pertur- 
bation theory: Such contributions, when included, alter the form of the static field 
action, but are not themselves statistically sampled. By this means it is possible 
to insure that higher-order approximations have a positive-definite weighting factor 
exp(— S e g), avoiding the sign problem inherent in a full quantum Monte Carlo. 

It is important to note that this approach has been explicitly designed to treat 
the low-temperature regime in which local moments exist, i.e., in which the most 
probable configurations of have non-zero amplitude (at least for some sites i). We 
do not address in this paper the problem of local moment formation per se. It is clear 
that at high temperatures the "free" £ 1 fluctuations have variance oc T, and 

it is not appropriate to use S e g(£) as a starting point for perturbation theory in this 
regime. Instead our purpose is to construct an approximation scheme which builds 
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in the presence of disordered local moments at temperatures satisfying 



tt 2 /U < T < U 



(37) 



in three dimensions and 



T < U 



(38) 



in two dimensions. 

It is also important to note that this approximation omits quantum effects as- 
sociated with "tunneling" of the order parameter fields, i.e., possible r-dependent 
instantons. As a result the Kondo effect is outside the present description. We re- 
strict our attention here to the implications of a classical (i.e., m = 0) order parameter 
description and do not discuss additional quantum effects. 

We next examine the behavior of S e s in the large-?/ limit, then comment on 
the action's relationship to other work in the literature. In the large-£7 limit the 
Hubbard model may be mapped to a quantum spin-| Heisenberg model with near- 
neighbor antiferromagnetic exchange integral J = 4t 2 /U. We show below that in 
this limit the approximate action S e g describes a classical Heisenberg model with 
exchange integral J = At 2 /U. This insures that in a two-dimensional system the 
approximation obeys the Mermin- Wagner Theorem, i.e., infinite-range spin ordering 
is absent down to T = 0. However this approach does not include the quantum 
fluctuations which renormalize the zero-temperature moment and reduce the finite- 
temperature correlation length. It is essential to include finite-frequency fluctuations 
in order to address these points. Nevertheless the static action S c fi achieves the goal 
of reproducing the exact classical critical behavior within a well-defined zeroth-order 
approximation. 

To demonstrate the mapping to the classical Heisenberg model, consider first the 
local one-site action (t = 0) in the large-f/ limit for (n) = 1: 



In this case there are only two eigenvalues e r . To find e r , choose the spin quantization 
axis along the direction of £. Then the field-dependent term in i?(£) is just 



-p^one— site 



m 2 - E l °s[ i + e ~^ (0 ] • 



(39) 



r 




(40) 



with eigenvalues 




a = ±1 . 



(41) 



Minimization of the action with respect to £ gives 




a 



(42) 



i.e., 



£ = \VlJ{m z ) , 



(43) 
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since 



(m \ = V a . (44) 

For large C/ the solution has 

(m z ) -> 1 

f - (45) 

As expected, the one-site gap is just 

ei-e^ = \U-{-\U) = U . (46) 

Note that this result has been obtained without breaking the rotational invariance 
of the one-site Hamiltonian, since the direction of £ is arbitrary, and the Boltzmann 
weight must be integrated over £. 

Thus to zeroth order in t/U, the half-filled system consists of fluctuating fields 
of fixed length VU/2 and uncorrelated directions. The first non-zero terms in 
perturbation theory enter at 0(t/U) 2 , and describe the process of one electron hop- 
ping to a neighboring site, then back to its starting point. Physically the hopping 
terms should induce correlations in the field directions at neighboring sites to lower 
the kinetic energy. Let the field directions at neighboring sites i and j be n % and W . 
Suppose also that n l may be obtained by rotating z in a right-handed sense by angle 
a about axis a _L z. Then the eigenstates with spin projection ±| along n % can be 
written 

\n\ T) = cos||2, t) + (-ia x + a y ) sin ^ \z, |) 

OL OL 

\n\ |) = (-ia x - a y ) sin - \z, t) + cos - \z, |) , (47) 



where 

The one-site solution gives 



+ <*t = 1 • (48) 



-v^7 = \\/U{n i ■ m' 1 ) , (49) 



so that I) will be occupied and |n\ |) will be empty (and higher in energy by 
U). Now apply second-order perturbation theory to the electron ground state for 
the frozen field configuration with direction vectors {n 1 }. The matrix element for 
hopping from 

\n\ T) |rP, T) - K I) |^, T) (50) 



is 



-t) (n j , I I n\ T) = (-*) (~i9 x - 6 y y sin | , (51) 
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where n j is obtained from n l by a right-handed rotation by 9 about axis 9 _L n\ The 
energy lowering of the ground state due to the virtual hop from i — > j and back is 
then 



-t)(-i9 x -9 y ) sin(fl/2) 
-U 



9 



U Sm 2 ' 



since 



+ 



(52) 



(53) 



The hopping from j i supplies the same energy, so the total perturbation energy 
lowering in this particular field configuration is 



-^E^ 2 (%/2) 

u (H) 



U 



]T(cos% - 1) 



(ij) 

u ^ 



(ij) 



u (ij) 



(54) 



— *. — # . 

where (S 1 )^ is the electron spin expectation value in configuration for t — 0. 

Note finally that the effective action at half-filling reduces to 



= /?Eiei 2 - log n [i + e '" £r( °] ■ 



(55) 



In the large-?/ limit 



exp[-^ (0] 



(56) 



with Eq the perturbed ground state energy, since all other occupancies result in an 
energy higher by at least U. Since /3J2i \C\ 2 i s a constant in this limit, the part of 
the action dependent on the fluctuating fields is precisely (3H eS , with 



(v) 



n ■ n J 



(57) 



establishing the mapping to the classical Heisenberg model. 

The behavior of the action away from half-filling in the large-?/ limit is not so 
easy to deduce, since the Hartree-Fock correction from the static charge interaction 
is crucial. If one assumes N e \ sites are singly occupied and A^ sites — iV el sites are 
empty, the one-site saddle-point equation implies = \/U /2 on the singly occupied 
sites and = on the empty sites. A degeneracy then arises, since the low- 
lying one-electron states on the singly occupied sites pick up a Hartree-Fock density 
contribution (U/2)(ri l ) = U/2, while the states on the empty sites have no Hartree- 
Fock correction. The two sorts of sites then have degenerate low-lying states in the 
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absence of t: The occupied state on each singly occupied site has energy — \U + \U = 

0, and the two states on each unoccupied site also have energy 0. As expected, the 
hopping becomes crucial for determining (ri 1 ) and the favored £ 1 configurations. Thus 
in the doped case there are important amplitude, as well as orientation, fluctuations 
in the S-H field configurations arising from S e s. 

We next comment briefly on the relationship of our approximation to other work 
in the literature. We have argued in Section 1 that perturbation theory in U fails 
to resum properly at low temperatures in the doped Mott insulator phase due to a 
fundamental change in the system's saddle-point structure. Nevertheless a number 
of authors 0-0 have observed that gaps or pseudogaps in the one-electron density 
of states for the positive-?/ and negative-?/ Hubbard models may be obtained within 
approximations which resemble perturbation theory about a non-interacting ground 
state. Such approximations have in common the feature that the one-electron self- 
energy S may be written in the form 

T 

E(k, iu n ) = — — 22 V(Q, itt m )G (k - Q, i(w n - fi m )) , (58) 

JVsites g m 

where V is a momentum- and frequency-dependent effective potential and Go is a 
bare electron propagator. Crucial points to note are that (i) the electron propagator is 
left bare; and (ii) higher order perturbation theory in V leads to a deterioration, rather 
than an improvement, in the approximation for the one-electron density of states. 
We argue that the reasons for this behavior are as follows: The local moments which 
underlie the Mott-Hubbard gap generate anomalous diagrams in perturbation theory, 

1. e., terms in the one-electron self-energy off-diagonal in momentum and spin. At 
zero temperature and half filling the moments have true long-range staggered order, 
so that a momentum-space description becomes possible. If the static off-diagonal 
self-energy is denoted Ax, with Q = (it, it), then the effective diagonal self-energy 
becomes 

iu n ) = A 2 ~G {k-Q) , (59) 
generating a Hartree-Fock antiferromagnetic gap of 2A~ in the density of states. In 



this case, as is well known, Go enters Eqn. |59] rather than G to avoid double count- 
ing. Thus approximations of the form in Eqn. |58| may be viewed as generalizations of 
Eqn. [5^ which "smear" the off-diagonal self-energy A over a range in momentum- 



and frequency-space |23|. In this sense these approximations are quite different from 
perturbation theory in V, and they do not provide support for the convergence of 
perturbation theory in the doped Mott insulator. We would argue that such approx- 



imations model only roughly the physics of this system p6fl . In the doped phase the 
order parameter field has strong non-Gaussian fluctuations in both amplitude and 
orientation, and these are coupled to the charge degrees of freedom. Furthermore, as 
discussed in reference |2l| , this general approach has difficulties addressing rare large 



amplitude fluctuations that play an essential role. 

We comment also on the relationship of the current approach to the phenomeno- 
logical spin-fermion (SF) model. In recent years a number of perturbative studies 
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27, 28, have been carried out for this model, which describes a system of elec- 
trons coupled to a Gaussian-distributed vector spin field S(Q, zO). The microscopic 
action S c s differs from Ssf in two crucial ways: (i) Ssf omits coupling to charge 
degrees of freedom, while S e s describes a highly nonlinear coupling between spin 
and charge. The form of this coupling is dictated by the Hubbard interaction itself, 
and its omission prevents a proper description of the doped Mott insulator phase, 
(ii) The finite-frequency corrections to the actions differ significantly, since we have 
introduced only a static S-H field to circumvent the problems associated with an 
equal-time decoupling. Finally, our point on the non-Gaussian nature of fluctuations 
in iSgff applies to Ssf as well, i.e., in the interesting parameter regime for this model we 
expect a conventional perturbative treatment to be impossible due to the formation 
of local moments. 



3 Evaluation of electron correlation functions 



Electron correlation functions for the action 5* in Eqn. have been calculated 
using the Metropolis algorithm for statistical sampling of £ * configurations. The Monte 
Carlo algorithm adopted relies on single-site updating and resembles a hybrid of 
Glauber and Kawasaki dynamics for the Ising or classical Heisenberg models. 
The algorithm is described more precisely in Section 4. 

Electron correlation functions have the same general form as in full quantum 
Monte Carlo calculations with continuous-valued S-H fields. Let A(c, c) be an elec- 
tron observable written in anticommuting c-number notation, and let 



(60) 



Then 



(A) 



J Dtj exp[-g (Q] J DcDc A(c, c) expj-S^, c, c)} 
J D£ DcDc exp[-So(0] exp c, c)} 



- I D£ exp[-S (0] 



/ DcDc A(c, c) exp[— Sf(£, c, c) 
/ DcDc exp[-£ f (£, c, c)] 



(61) 



where 



Z f (£) = J DcDc exp[-S f (£, c, c)] 



nti+ 



,-/M0i 



(62) 
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and the last line in Eqn. 61 follows by multiplying and dividing by Zi within the 
integral. The expression in square brackets is just an electron thermal average in the 
presence of the background £ 1 configuration. Thus finally 



1 



Di exp[-5 eff (0] (A{c, c)) ( 



where 



((Mc, c)> € > 



5 (0 - logZ f (0 



(63) 



(64) 



The simple form for S e s means that approximations which incorporate finite-frequency 
fluctuations, i.e., which generalize Zf(£), may in principle be developed to extend the 
approach worked out in detail here (see Section 6). 

We derive below the form of several one- and two-body electron correlation func- 
tions (A(c, c))^ in the presence of the background fields. As stated previously, 



= E G aa (ii; t = ) $ 

a 

= J2f(er)\(r\ia)\\ 



where 



G aa >(ij; r) € 



(c fa (r)^(0)) < 



(65) 



(66) 



in terms of anticommuting c-numbers; and the |r) are eigenstates of the spin-dependent 
Hamiltonian H(£) in Eqn. Likewise, the spin polarization at site i is 

aa' 

= E ° aa ' G a 'a(ii; T = 0~)| 
aa' 

= E /( £ r) ^ (r\ia) (ia'\r) , (67) 



in terms of the S^-basis states \ia). Writing out the vector components explicitly 
gives 



(m*) f 



E/(^)[(U|r)(r|U) + (^t|r)(r|z|) 

r 

E /(^) i l r ) ( r N t) + «(« T \r) (r\i |) 

r 

E/c^fi^tior - i(uir)i 2 



(6? 
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Note that these last results generalize to give the one-body propagator 
GiajpiiuJt = [ drd^ [_( c -( r )^(0))/ 

= ^ (ia\r) (r\j(3) (gg) 

The real-axis spectral density may be calculated directly by analytic continuation 
before the field averaging is performed. Thus the one-particle spectral density for 
spin a at site i after field-averaging is 

N ia {uj) = -- (lmG aa (ii; u + iO + )t: 

7T x 

^\(ia\r)\ 2 5(u;-e r )). (70) 

r 

The total one-particle spectral density for spin a is 

N a (uj) = E N ia (uj) . (71) 

i 

(Note that in a properly equilibrated calculation translational and spin-rotational 
invariance are restored by field-averaging, and the quantity iVj a is site- and spin- 
independent.) The momentum-resolved spectral density Ak a {uj) is only slightly more 
difficult to calculate. Note again that translational invariance is only restored by 
field-averaging, so that the propagator in a single field configuration is not diagonal 
in momentum. After averaging we have 

A ka (^) = ( E % - e r ) • t^- E e-***-^ (ia\r) (r\ja) ) (72) 

r -''sites jj 

All the spectral densities may be calculated conveniently by histogram binning, 
i.e., by incrementing a bin centered on frequency u>i whenever an eigenvalue e r falls 
in the range [cui — (l/2)Au, u>i + (l/2)Au). In this way spectral densities may be 
obtained which integrate exactly to unity within machine precision. 

In some cases it may be desirable to compare imaginary-time propagators Gk a { T ) 
with the results from a full quantum Monte Carlo simulation. The necessary expres- 
sions within the present approach are 



. " Er e- £rT [1 - /(e r )l F r ka ), r > 



where 



Ka = T7 — E e-*<*-™ (ia\r) (r\ja) . (74) 

-< * sites 
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A quantity of particular interest in the local moment regime is the average double 
occupancy 

(D) = T^-E(We). (75) 

-''sites j 

where 

(D% = (nV^ . (76) 

In a general £ 1 configuration with broken spin-rotational invariance the appropriate 
expression for (D 1 )^ is not (n^)g (ft )<•■ Instead 

(£>*)* = lim (c iT (r)c iT (r)c i - L (^ , )c ii (^))€ 

= (n iT ) ? (n il )t + (c^c il )^ (&c% . (77) 

Note 

<c*V^ = l[K) e + i(m^] 

(c l V^ = -i[K.) 5 -z(m;) ? ], (78) 

while 

§[<»*>« + K)d 

- K)d • (79) 
Thus 

(£>*> € = \[{n% - . (80) 

The double occupancy assumes an explicitly spin-rotation-invariant form as required. 

The expression in Eqn. |5l] should in principle be supplemented by zero- and 
finite-frequency RPA corrections to insure that Eqn. [77| reproduces the r — > r' limit 
of the general two-body correlation function. (As noted below in the discussion of 
response functions, RPA corrections appear at zero frequency because the charge 
interaction is being treated within Hartree-Fock. There is no zero-frequency RPA in 
the spin channel, since that portion of the interaction is being treated exactly.) RPA 
corrections in equal-time expectation values (n ia {T)n %a (r)) are actually problematic 
in any perturbation theory for the following reason: The Pauli Principle implies an 
exact identity 

(n ia (r)n ia (r)) = (n ia ) . (81) 

In perturbation theory, this identity can only be preserved exactly in approximations 
with a crossing-symmetric two-body vertex. Such a vertex appears in order-by-order 
perturbation theory (and in parquet theory), but not in a standard RPA, whether 
about the unperturbed vacuum or about broken-symmetry states as in the present 
case. The RPA can be augmented to make it crossing-symmetric by incorporating 
the missing contributions in the crossed channel (see Figure 3), but the result is no 
longer conserving in the Baym-Kadanoff sense. The RPA corrections to the double 



(n% = 
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occupancy are in any event expected to be small in the present case. Since they are 
also problematic, we omit them from the calculation of (D) reported here. 

We next consider the static spin and charge susceptibilities Xs(Q) an d Xc(Q)- As 
mentioned above, one expects the appearance of RPA corrections in the static charge 
response, since the charge interaction is being treated by "slave" Hartree-Fock, but 
naively one might expect no corrections in the spin response. This is not correct: An 
average must be performed over the symmetry-breaking field £ l , and in the presence 
of spin and charge excitations are mixed. For this reason the repulsive charge 
interaction makes its presence felt weakly in the spin channel, as well as in the charge 
channel. 

The static spin response function takes the form 



(82) 



X% = d((m^)/dhi 

h=0 

where the assumed local coupling to the external field h? is 

AH = -W • £ c] a a aa ,c ja , . (83) 

aa' 

Temporarily neglecting RPA corrections, the field derivative yields 

X% = ~ ( H> € ) ( (mi), ) + C<X«' E ) . (84) 

1 1 n n > ? 

The first term on the right-hand-side vanishes for zero applied field (since the system 
has no spontaneous magnetization at finite temperature), and the second term (see 
Figure 4) reduces to 

Xj v = ^(( m ^)(( m l)() + (~ T E 2 G pa{ji\ i^n)^G a ,p,{ij; iUnhVaa'Vp'p) , 

(85) 

i.e., formally an "anomalous" and a "normal" contribution. The anomalous contri- 
bution may be evaluated using previous expressions for the site spin polarizations, 
while the normal term may be written as the configurational average of 

where 

M^ iaM (i^ = -EO'Ak) (r|< ai > (<a 2 | a > (s\j(h) f *f r l~ f ^ ■ (87) 

(Note that for e r = e s , the ratio involving the Fermi functions must be replaced by 
/' ( £ r)<W ) 
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The static charge response function is 
where in this case the assumed coupling to the external field is 

AH = ~fi J 4* C ja ■ 



39) 



In complete analogy with the spin response calculation the result before the inclusion 
of RPA corrections is 

1 ~ ' (n% (ni)t ) - (n) s 



T 



+ E 

aa'(3(3' 



(90) 



with the normal "bubble" M defined as in Eqn. 

As emphasized above, RPA corrections necessarily appear when an interaction is 
treated perturbatively by Hartree-Fock. In the present case, where spin and charge 
excitations mix before the configurational average is performed, RPA charge vertex 
corrections appear in both the static charge and spin response functions. 

The calculation of RPA corrections for Unite-frequency response functions is trick- 
ier than in conventional perturbation theory because the zero-frequency spin ver- 
tices have been explicitly removed from the theory. This affects response functions 
at iQ 7^ by partially removing vertices which carry zero frequency in a crossed 
particle-hole channel. The analytic continuation to real energies is also complicated 
by this feature of the theory. We defer discussion of spin and charge response at fi- 
nite frequency and concentrate in the present paper on RPA corrections to the static 
susceptibilities. 

The repulsive RPA charge vertex produces a geometric series of corrections as in 
conventional perturbation theory. Due to the lack of symmetry in individual config- 
urations, the series can only be summed in matrix form. Writing 



y 0102, /3i/32 



(91) 



2 U 5ij5 aia2 5 i3 2 i3\ 1 

the complete RPA-corrected expressions for the spin and charge susceptibilities are 
given by 



Y 13 



and 



1 

T 



{ml) 1 



( [Mil + VMy 1 



aa'/3f3' 



aa',(3f3> 



{iQ = o) 5 



[A v 
a aa' <J f3'0 



(92) 



X 



1 

T 



+ E 

aa'/3/3' 



~ (n) 2 

[m(i + vm)- 1 



act', (3(3' 



0) c )<W<W • (93) 
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Finally the momentum-resolved spin susceptibility takes the form 

x,AQ) = t^-E^-^x^, (94) 

-'''sites jj 

and similarly for the charge susceptibility. (Measurements of the RPA-corrected 
susceptibilities require the inversion of small matrices, but the effect on overall cal- 
culational efficiency is negligible.) 



4 Description of the Monte Carlo algorithm 

The effective action in Eqn. [64] has been simulated using the Metropolis algorithm. 
In comparison with a classical Heisenberg spin model, the present theory has several 
complications: (i) The theory has three degrees of freedom (orientation on the unit 
sphere and amplitude) at each site, rather than two. (ii) Local field updates induce 
global changes in the action through the Fermi partition function Zf(£). (iii) A 
Hartree-Fock convergence step is necessary after each spin field update to keep the 
"slave" charge fields (n 1 )^ on the saddle point, (iv) Single-site updating (i.e., Glauber 
dynamics) is insufficient to achieve efficient equilibration for (n) ^ 1. An analog of 
Kawasaki dynamics, or spin swapping, is required, since the doping 1 — (n) plays a 
role similar to a conserved order parameter. 

We consider each complication in turn, (i) The natural physical interpretation 
of field fluctuations in terms of amplitude and orientation suggests an optimal rep- 
resentation for the functional integral over field configurations. At half-filling it is 
expected that the only important fluctuations for T —>■ correspond to magnons, i.e., 
slow modulations of the field orientation for fixed amplitude. A Cartesian represen- 
tation of the fields (£* , £*, £*) with independent component updates leads to a small 
acceptance ratio in the low-temperature limit. Therefore we write 



/ 



•2tt 

dCdCdC = I (r l )' 2 dr l I d(cos9 l ) I d<j>\ (95) 



choosing as independent variables r l , x l = cosO 1 , and <ft l . Variations in the measure 
(r*) 2 at each site are incorporated into the Metropolis updating factor, and moves 
are restricted, in a way consistent with detailed balance, to lead to r % > 0. Proposed 
changes to r % are restricted to a maximum scale Ar, while large-scale changes are 
permitted in x l and <p l . 

(ii) For the initial Monte Carlo calculations reported here we have simply diago- 
nalized the effective Hamiltonian matrix H{£), which has row dimension 2iV S it es , for 
each proposed £ configuration. The time for a single site update thus grows as the 
cube of the system size, and the time for a lattice sweep as the fourth power of the 
system size. We believe that a hybrid molecular-dynamics/Monte-Carlo algorithm 
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can be developed to significantly reduce the time needed for the simulation of large 
lattices J31] . 

(iii) The Hartree-Fock treatment of the static charge interaction is an integral com- 
ponent of the present approach. More sophisticated treatments of fluctuations would 
also in principle involve a self-consistent calculation in the presence of symmetry- 
breaking background fields. For the parameter sets studied in the present paper, the 
Hartree-Fock iteration loop generally converges in three to four steps. Since each step 
requires a matrix rediagonalization, overall calculational time scales with the average 
number of Hartree-Fock iterations. 

(iv) At half-filling the Hartree-Fock charge-field convergence becomes trivial, 
and (n l ) remains equal to unity at all sites in all configurations. This just reflects 
the particle-hole symmetry of the model for this special parameter set. Away from 
half-filling the situation is quite different. Typical configurations contain a few sites 
with small fields, while the majority of sites have fields with amplitudes close to the 
half-filled value. Conventional single-spin updating tends to leave "holes" trapped 
at isolated locations, since moving a hole requires passing through high-energy in- 
termediate configurations. In this sense the number of hole 1 — (n) functions like 
a conserved order parameter, preventing the equilibration of the system by Glauber 
dynamics. To combat this effect we have introduced a second type of updating move, 
viz., a swap of the spin fields £ at two sites chosen at random. When the spin fields 
are swapped, the "slave" charge fields are dragged along, and the system equilibrates 
to a uniform set of site occupancies ( (ri 1 )^), even though the hole distribution in 
individual configurations is highly non-uniform. 



5 Numerical results for a small Hubbard lattice 

The results for 4x4 periodic lattices reported below were obtained on an IBM 
PowerPC work station and a Sun-3000 server. For all runs reported here, a total of 
100 Monte Carlo warmup sweeps and 5000 sampling sweeps were carried out. Runs 
on 8 x 8 and 10 x 10 lattices are now in progress on a Sun-10000 work station array 
at the San Diego Supercomputer Center. 

The total density of states for spin a, N a (w), is plotted in Figure 5 for U/t — 4 and 
in Figure 6 for U/t = 8, both at temperature T/t = 0.125. The chemical potential 
n/t is lowered from U/2 (half-filling) as holes are doped into the insulating state. 
The mean site occupancy (n) and double occupancy (D) are summarized in Table 1. 
Note that the double occupancy at half-filling is reduced by factors of order two and 
six from the value 0.25 for a non-interacting system. A prominent Mott-Hubbard 
gap, or pseudogap, appears in the density of states for both parameter sets at half- 
filling (Figures 5a and 6a). The remnants of the discrete quantum energy levels in 
the non-interacting system remain visible. Note again that these systems are not 
magnetically ordered: The spin-up and spin-down densities of states are identical 
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within statistical sampling errors, and the spin susceptibility is isotropic, as discussed 
below. 

As the chemical potential is pulled below /i = U/2, the site occupancy changes 
slowly at first, particularly for the large charge gap which appears when U/t = 8 
(Figures 6b-c). This is just what one would expect for rigid-band doping in a semi- 



conductor, and such behavior has been observed in previous Lanczos |32L |33| and 



quantum Monte Carlo studies of the Hubbard model on small lattices [35|| . Note 
that the rigid-band picture begins to break down as soon as the site occupancy 
drops appreciably below unity (Figure 5b). In the doped system with U/t — 4 and 
fi/t = 1.2, spectral weight has been drawn from the lower edge of the upper Hubbard 
band into the gap. Due to this transfer of weight, the gap feature rapidly collapses 
with doping, though remnants remain for a wide range of fillings (Figures 5c-d). By 
the time \i/t reaches 0.4 for U/t = 4, the shape of the non-interacting band has begun 
to re-emerge, though levels are still significantly broadened by the fluctuating spin 
background. 

This broadening is especially evident in the evolution of the k-resolved spectral 
density for points near the Fermi surface. Recall that a special symmetry of the 
4x4 lattice induces an equivalence between the generally distinct k-points (tt, 0) 
and (tt/2, tt/2). This means all points on the half-filled Fermi surface have the same 
Ak a (uj) for the small systems discussed here. To contrast the behavior of points near 
and far from the Fermi surface, we plot Ak a (u) for k = (tt, 0) and k = (tt, tt) at 
U/t = 4 in Figures 7 and 8. While the (tt, tt) point at the top of the band (Figure 
8) has a sharp quasiparticle-like spectral density whose profile changes only slightly 
with doping, the (tt, 0) point behaves quite differently. At half-filling Ak a for this 
point is symmetrically split between the upper and lower Hubbard band edges, just 
as one would expect in an antiferromagnetic Hartree-Fock solution |T]]. Away from 
half-filling, the spectral density remains exceptionally broad, reflecting the continued 
presence of strong short-range antiferromagnetic correlations between the fluctuating 
spins in the doped system. As noted above in the discussion of N(u), remnants of the 
half-filled gap remain for a wide range of fillings. The spectral broadening decreases 
for fillings far from unity, consistent with the restoration of the non-interacting band 
profile (and presumably the recovery of Fermi liquid behavior). 

A broadening of the spectral density is also present in perturbative studies, includ- 
ing self-consistent-field methods such as FLEX || ; however, the details of the broad- 
ening mechanism are quite different in the present case. FLEX and phenomenologi- 
cal spin-fluctuation models [p9j j36fl assume a small-amplitude scattering mechanism 
which may be treated by perturbation theory about a non-interacting (i.e., Fermi 
liquid) background. In contrast the present approach treats scattering by large- 
amplitude background distortions, with highly non-Gaussian static fluctuations in 
both amplitude and orientation; these distortions evolve smoothly into the ordered 
antiferromagnetic background at half filling and zero temperature. 

The behavior of the spectral density at the (tt, 0) point for U/t = 8 is illustrated 
in Figure 9. As expected, the spectral weight is split symmetrically at half-filling, and 
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the two peaks are separated by an energy of order U. In contrast with the U/t = 4 
result, the weight is distributed broadly throughout the upper and lower Hubbard 
bands. Signs of a more narrow coherent feature in the doped lower band are present 
for fi/t = 1.8 (Figure 9c). 

The momentum dependence of the static charge and spin susceptibilities is plotted 
in Figures 10 and 11 for U/t = 4 and 8 at temperature T/t = 0.125. Note that the 
uniform charge susceptibility (or compressibility) is exponentially small at half-filling, 
as expected for a system with a charge gap. In addition the compressibility increases 
very rapidly with doping: Note particularly the increase for U/t — 8 (Figure 11a). For 
general fillings Xc{Q) is relatively small (note the difference in scales for the charge and 
spin susceptibility) and structureless at this temperature. This is expected in such a 
small system. Possible manifestations of charge ordering within the fluctuating spin 
background may become apparent, however, in larger systems. 

In contrast the spin susceptibility is large and strongly peaked at the antiferromag- 
netic wave vector (n, ir). The peak is expected to move off (ir, ir) for doped systems 
on larger lattices. For the parameter sets illustrated here, the spin susceptibility is 
dominated by the contributions of the "anomalous" term 



The [iv = zz component of Xs is actually plotted in Figures 10 and 11. For 5000 
lattice sweeps the agreement between the xx, yy, and zz components of the (ir, ir) 
spin susceptibility is at the level of the statistical error bars, i.e., a few percent. (The 
agreement can be arbitrarily improved by increasing the Monte Carlo sampling time.) 



6 Extension to higher order Baym— Kadanoff approximations 

The present approach consists of simulation of the action for a set of local static 
fields coupled to electron spin degrees of freedom. When the action is constrained to 
be at its temperature-dependent global saddle point (if such a point indeed exists in 
the most general case), these static fields are determined by inhomogeneous Hartree- 
Fock equations and reduce to the anomalous fields of conventional diagrammatic per- 
turbation theory. Just as Hartree-Fock theory with an anomalous self-energy may be 
extended using the Baym-Kadanoff prescription [P^| , there is no formal barrier pre- 
venting such an extension in the present case. The computational cost of carrying out 
self-consistent -field (SCF) perturbation theory in a broken-symmetry background is 
quite large, however, since total momentum and spin are no longer good quantum 
numbers. In this section we sketch out the steps necessary for a generalized approach 
and note several of the technical pitfalls which must be studied and overcome. 




(96) 
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The form for the SCF action must be determined separately for each configuration 
of the fields. In full analogy with Baym-Kadanoff theory for a uniform system, 
the action takes the form 

S SCF (0 = /?£ \C\ 2 - TV log (-G^T 1 + fiT, (97) 

i 

where 



and 



= -Tr(SG b ^) + $ (98) 
5$ 



In these formal expressions G^ CF is the discrete-time propagator (retained to insure 
a proper high-frequency regularization of the Tr log term); 

x = (r*, r) , (100) 

and $ is the SCF generating functional (e.g., a sum of ring and ladder diagrams in 
FLEX). 

Since the SCF action may be interpreted as /3JF SCF , with JF SCF a well-behaved 
free energy, the statistical weighting factor exp(— 5' SCF ) is positive definite. This is in 
sharp contrast to the behavior in imaginary-time quantum Monte Carlo simulations, 
where the weighting factor can be positive or negative with nearly equal probability. 

What is not so clear is the nature of approximations which insure the stability of 
S , particularly for exponentially rare field configurations. The interaction which 
enters S SCF is not the original Hubbard interaction, since the zero-frequency spin 
components have been removed for special treatment. This means that there are 
no obvious zero-frequency instabilities left to appear in S SCF . Behavior at finite 
frequencies must also be addressed, however, and we are not prepared to comment 
on the stability of general Baym-Kadanoff approximations (and in particular FLEX) 
in the present paper. 

As an example of the Baym-Kadanoff prescription in practice, note that when the 
static charge interaction is treated within Hartree-Fock 

E aal (xx') = \U(n l )^8 iV 8 TT >8 a(Tl . (101) 

The appropriate generating functional is 

* = J^E <»*>€■ (102) 

i 

Furthermore 

Tr(SG SCF ) = 2$ (103) 

and so 

7 = -$ 

= -^E <»*>!■ (104) 
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Note that this form for the Hartree-Fock "constant" in the action is exactly that 
which has appeared previously in Eqn. [H]. 

7 Some speculations and conclusions 

We believe the present Hubbard model analysis, based on Monte Carlo simula- 
tion of a classically disordered local spin field, incorporates several desirable features 
which have eluded previous perturbative studies. Most important the present theory 
has classical critical behavior at half-filling and describes the formation of a Mott- 
Hubbard gap (or pseudogap) in a system with a finite spin correlation length. In 
addition the spectral densities for the doped systems have the features suggested by 
exact methods. 

The present static field analysis can be extended to significantly larger systems. 
Since the computations involve Monte Carlo field sampling and have reasonably small 
memory requirements, it would be natural to carry out studies on multiple processors 
with different random seeds, then to merge the equilibrated results from numerous 
short runs. 

It is almost unavoidable to speculate on the implications for charge density wave 
or stripe formation |37|, |38| and d-wave singlet pairing in the cuprate superconductors. 
Since the present approach reduces to standard inhomogeneous Hartree-Fock for T — > 
0, Hartree-Fock results previously obtained for stripe formation should re-emerge 0]. 
By allowing sampling of configurations removed from the zero-temperature saddle 
point in a well-defined way, the present approach provides a framework to describe 
the destruction of striped order by thermal fluctuations. 

To study potential d-wave superconductivity it seems essential to incorporate the 
exchange of low-frequency spin fluctuations about the disordered background con- 
figurations. The low-energy fluctuations which become magnons in the directions 
transverse to an ordered spin state (i.e., the Goldstone modes of the order param- 
eter field) must evolve into a broadened spectrum of fluctuations in the important 
disordered configurations. If low-frequency spin fluctuations are important, the pair 
order parameter should exhibit retardation effects, in contrast with the spin order 
parameter. In addition the disordered background has important implications for 
the onset of pairing. Presumably the pair field appears initially with a highly non- 
uniform distribution in a small subset of £ configurations. Such behavior would imply 
superconducting precursor effects in the density of states well above the transition 
temperature for global phase coherence. The effects would be more akin to the behav- 
ior of a system of superconducting grains with variable grain size than to conventional 
Kosterlitz-Thouless phase fluctuations in a clean two-dimensional superconductor. 

In any case the scenario described above is at this stage purely speculative. The 
next steps in the development of the present approach consist of (i) studies of larger 
systems; (ii) analysis of the spin and charge response at finite frequencies, with partic- 
ular care given to the treatment of potential instabilities in rare field configurations 
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(see Section 6); and (iii) analysis of pairing eigenvalues and eigenvectors from ex- 
change of fluctuations in sets of non-uniform configurations. 
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Table 1. Mean occupancy (n) and double occupancy (D) for parameter sets at T/t = 0.125. Statis- 
tical errors are in the last digit. 
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Figure 1 : Diagrammatic representation of the Hubbard interaction vertex. 
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Figure 2: Rewriting of the Hubbard interaction vertex in the two crossed particle-hole channels, 
(a) Transverse spin channel, (b) Longitudinal spin and charge channel. 
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Figure 3: Example of reducible and irreducible particle-hole vertex diagrams related by crossing 
symmetry. The third-order reducible vertex appears in an RPA summation. It must be augmented 
by the irreducible particle-hole vertex shown in order to preserve the Pauli Principle identity in 
Eqn. 
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Figure 4: Contributions to the static spin response function, (a) "Anomalous" diagram, (b) "Nor- 
mal" diagram. 
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Figure 5: Total density of states for spin a, N a (uj). The interaction strength is U/t = 4, and the 
temperature is fixed at T/t = 0.125. Results are shown for four different chemical potentials. See 
Table 1 for the mean site occupancy and double occupancy values. Note that in the absence of 
statistical sampling error Nf(ui) = N±(u>). Therefore the measured difference in these two quantities 
provides a point-by-point error estimate. The vertical dashed line at cj = is drawn to emphasize 
the position of the chemical potential, (a) fi/t = 2.0. (b) fi/t = 1.2. (c) fi/t = 0.8. (d) fj,/t = 0.4. 
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Figure 6: Total density of states for spin a, N a (ui). The interaction strength is U/t = 8, and the 
temperature is fixed at T/t = 0.125, as in Figure 5. See Table 1 for the mean site occupancy and 
double occupancy values, (a) [ijt = 4.0. (b) [ijt = 2.4. (c) [ijt = 1.8. 
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Figure 7: Momentum resolved spectral density Ak a (uj). The model parameters arc U/t = 4 and 
T/t = 0.125, as in Figure 5. The momentum point chosen is k — (tt, 0). (a) /i/i = 2.0. (b) jj,/t = 0.8. 
(c) fi/t = 0.4. 
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Figure 8: Momentum-resolved spectral density Ak a (^)- The model parameters are as in Figure 5 
and 7. The momentum point chosen is k = (w, n). (a) fi/t = 2.0. (b) fi/t = 0.4. 
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Figure 9: Momentum resolved spectral density Ak a (uj). The model parameters arc U/t = 8 and 
T/t = 0.125, as in Figure 6. The momentum point chosen is k — (ir, 0). (a) [ijt = 4.0. (b) fi/t = 2.4. 
(c) fi/t = 1.8. 
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Figure 10: Static charge susceptibility Xc(Q) an d spin susceptibility Xs(Q) f° r U/t = A and T/t = 
0.125. The susceptibilities are plotted along the triangular Brillouin zone contour from r — > X — > 
M — > T, i.e., from (0, 0) — > (7r, 0) — > (7r, 7r) — > (0, 0). Results are shown for the chemical potentials 
and site occupancies employed in Figure 5. Unless shown, statistical error bars are smaller than the 
plotting symbols, (a) xAQ)- (b) Xs(Q)- 
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Figure 11: Static charge susceptibility Xc(Q) and spin susceptibility Xs(Q) for U/t — 8 and T/t = 
0.125. The susceptibilities are plotted along the triangular contour used in Figure 10. Results are 
shown for the chemical potentials and site occupancies employed in Figure 6. (a) Xc(Q)- (b) Xs(Q)- 
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